Baseline exploration

p <- df %>% 
  select(run_number = .run.number., chance_imp = importance.of.chance, 
         step = .step., contains("gini")) %>% 
  pivot_longer(contains("gini")) %>% 
  filter(factor(chance_imp) %in% factor(seq(0, 1, by = .1)),
         !is.na(value)) %>% 
  mutate(name = str_remove(name, "gini_") %>% str_to_title()) %>% 
  ggplot(aes(step, value, colour = factor(chance_imp))) +
  geom_smooth() +
  facet_wrap(vars(name)) +
  scale_colour_viridis_d(option = "C", alpha = .5,begin = .1, end = .9) +
  labs(y = "Gini index", 
       colour = expression(paste("Randomness parameter ", italic("c"))),
       caption = "n = 100 agents; summarised from 20 runs") 
p

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

plotly::ggplotly(p)

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

Plot of detail

p <- df_detail %>% 
  select(run_number = .run.number., chance_imp = importance.of.chance, 
         step = .step., contains("gini")) %>% 
  pivot_longer(contains("gini")) %>% 
  mutate(name = str_remove(name, "gini_") %>% str_to_title()) %>% 
  filter(!is.na(value)) %>% 
  # mutate(chance_imp = factor(chance_imp, labels = chance_imp %>% unique() %>%
  #                              as.numeric() %>% scales::percent())) %>%
  ggplot(aes(step, value, colour = factor(chance_imp))) +
  geom_smooth() +
  facet_wrap(vars(name)) +
  scale_colour_viridis_d(option = "C", alpha = .5, begin = .1, end = .9) +
  labs(y = "Gini index", 
       colour = expression(paste("Randomness parameter ", italic("c"))),
       caption = "n = 100 agents; summarised from 20 runs") 
p

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

Distribution of publications

df_end_selected <- df_end %>% 
  select(run, importance.of.chance, gini_publications, n_publications,
         total_grants)

p_density <- df_end_selected %>% 
  ggplot(aes(n_publications, fill = factor(importance.of.chance), group = importance.of.chance)) +
  geom_density(alpha = .5, adjust = .6) +
  scale_fill_viridis_d(option = "C") +
  geom_curve(
    aes(x = 2500, y = .005, xend = 3000, yend = .001), curvature = -.35,
    arrow = arrow(length = unit(0.03, "npc"))
  ) +
  annotate("text", x = 2000, y = .005, 
           label = expression(paste("Bimodal distribution\nfor low values of", 
                                    italic("c")))) +
  labs(x = "# of publications per research group",
       fill = expression(paste("Randomness parameter ", italic("c"))),
       y = "Density") +
  theme(legend.position = "top") +
  guides(fill = guide_legend(byrow = TRUE)) 
p_density

ginis <- df_end_selected %>% 
  distinct(importance.of.chance, gini_publications) 

p_ginis <- ginis %>% 
  ggplot(aes(gini_publications, factor(importance.of.chance),
             fill = factor(importance.of.chance))) +
  geom_boxplot(width = .5, show.legend = FALSE, alpha = .5) +
  geom_jitter(height = .1, width = 0, show.legend = FALSE, alpha = .5) +
  scale_fill_viridis_d(option = "C") +
  labs(x = "Gini index of publication-distribution",
       y = expression(paste("Randomness parameter ", italic("c")))) 

p_density + p_ginis +
  plot_layout(heights = c(3, 2)) +
  plot_annotation(tag_levels = "A")

df_end_selected %>% 
  filter(n_publications < 1500) %>% 
  ggplot(aes(n_publications, fill = factor(importance.of.chance), group = importance.of.chance)) +
  geom_density(alpha = .5, adjust = 1) +
  scale_fill_viridis_d(option = "C") +
  labs(x = "# of publications per research group",
       fill = expression(paste("Randomness parameter ", italic("c"))),
       y = "Density") +
  theme(legend.position = "top") +
  guides(fill = guide_legend(byrow = TRUE)) 

df_end_selected %>% 
  ggplot(aes(n_publications, fill = factor(importance.of.chance),
             group = importance.of.chance)) +
  geom_histogram()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

df_end_selected %>% 
  ggplot(aes(x = n_publications, y = factor(importance.of.chance), group = importance.of.chance)) +
  geom_density_ridges(scale = 2)

## Picking joint bandwidth of 43.2

p <- df_end_selected %>% 
  filter(n_publications < 1500) %>% 
  ggplot(aes(x = n_publications, y = factor(importance.of.chance), 
             group = importance.of.chance, fill = stat(x))) +
  geom_density_ridges_gradient(scale = 1.7, show.legend = FALSE) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_discrete(expand = c(0, 0)) +
  scale_fill_viridis_c(option = "C", direction = -1, end = .9) +
  coord_cartesian(clip = "off") + 
  hrbrthemes::theme_ipsum_rc(base_family = "Hind",
                             grid = FALSE) +
  # theme_ridges(grid = FALSE, center_axis_labels = TRUE, font_family = "Hind") +
  labs(x = "# of publications per research group",
       y = expression(paste("Randomness parameter ", italic("c"))),
       caption = "x-axis is truncated at 1500")
p

## Picking joint bandwidth of 25.4